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1. Introduction 

Extrasolar planetary systems have became a major challenge for contemporary astrophysics 
and dynamical astronomy. One of the most difficult problems in this field concerns the orbital 
stability of such systems, in particular when related to the observations and their interpretation. 
Usually, the investigations of long-term evolution are the domain of direct, numerical integra- 
tions. The stability of planetary systems is often understood in terms of the Lagrange definition 
implying that orbits remain well bounded over infinite time. Other definitions may be formu- 
lated as well, like the astronomical stability (Lissauer||1999|) requiring that the system persists 



over a very long, Gyr time-scale or Hill stability ( Szebehely||1984| ) that requires the constant 
ordering of the planets. In our studies, we prefer a more formal and stringent definition related 
to the fundamental Kolmogorov-Amold-Moser Theorem (KAM), see Arnold ( 1978 ). Planetary 
A^-body systems, involving a dominant mass of the parent star and significantly smaller planetary 
masses, are well modeled by close-to-integrable, Hamiltonian dynamical systems. According to 
the KAM Theorem, their evolution may be quasi-periodic (with a discrete number of fundamen- 
tal frequencies, which are stable forever — and also stable with respect to the other notions of 
stability quoted above), periodic (or resonant; stable or unstable) or chaotic (with a continuous 
spectrum of frequencies, and unstable). In the last case, initially close phase trajectories diverge 
exponentially, i.e., their Maximum Lyapunov Characteristic Exponent (MLCE, denoted also by 
cr) is positive. This understanding of the global structure of the phase space is widely adopted, 
in particular with regard to the Solar system dynamics (e.g.. W isdom & Holman^^ 1991; Laskar 



1990 



2001 



2002 



Holman & Murray! 1996||Malhotra|l 9981 [Nesvorn y & M orbidelli|1999 



Robutel & Laskar 



Murray & Hol manf2001t[Michtchenko & Ferraz-Mello|200H|Lecar et al.|2001[ 



Morbidelli 



and references therein.). However, the distinction between regular and chaotic trajectories 
is a difficult task that, in practice, may be resolved only with numerical methods relying on 
efficient and accurate integrators of the equations of motion. 



2. Stability indicators 

To detect chaotic motions in the phase space, many numerical tools are available. Concerning 
the dynamics of close to integrable Hamiltonian systems, these tools can be roughly divided onto 
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two c lasses: spectral algorithms that resolve the fundamental frequencies and/or their diff usion 
rates (|Laskar||l993l ISidhchovsky & Nesvorny||l997l iMichtchenko & Ferraz-Mello||200l|), and 



methods based on the divergence rate of initially close phase trajectories, expressed in terms of 
the Lyapunov exponents (Benet tiiTet al.|1980{|Froeschle|1984| ). These fast indicators can be cor- 
related with geometrical evolution of orbital osculating elements, like the maximal eccentricity 
(maxe), maximal amplitude of the critical angle of a resonance (max6>) or with event time Te 
(indicating a collision or an ejection of a body from the system). Unfortunately, no general rela- 
tion between these indicators can be determined ( [Lecar et al.|2Q01t[Michtchenko & Ferraz-Mello| 



2QQ1| ). The event time Te, relying on CPU-intensive long-term integrations of the equations of 



motion ( [Lecar et al.|2001| ) can be considered as a direct measure of the astronomical stability. 



In our work, among the the spectral tools, we often choose the method invented by |Michtchenk4 
& Ferraz-Mello| ( [200T] ); its idea is very simple — to detect chaotic behavior one counts the num- 



ber of frequencies in the FFT-spectrum of an appropriately chosen dynamical signal. We deal 
with conservative Hamiltonian systems; so in a regular case, the spectrum of fundamental fre- 
quencies is discrete and we obtain only a few dominant peaks in the FFT spectrum. Chaotic 
signals do not have well defined frequencies, and their FFT spectrum is very complex. The num- 
ber of peaks in the spectrum above some noise level p (typically, p is set to a few percent of the 
dominant amplitude) tell us about the character (regular vs chaotic) of a given phase trajectory. 

The basic tool to discover exponentially unstable bounded orbits, i.e. chaotic orbits, is the 
Maximum Lyapunov Characteristic Exponent (MLCE) cr. The direct computation of the MLCE 
is based on the analysis of the tangent vectors 8 which are solutions to the variational equations 
of the equations of planetary motions: 

= fix), 

where x denotes the state vector including coordinates and momenta, and / stands for the grav- 
itational forces. For its solution = (f){t) we define: 

^^S = A{t)S, A{t):=^m)], S=\\S\\. (2.1) 
Asymptotically, the MLCE value is given by j Cincotta & Sim6|2QQQ| ): 

CT = lim i ^d8. (2.2) 

If a converges to some positive value, we conclude that the nominal orbit (f) and some initially 
close orbit diverge exponentially at the rate exp(crt). Two practical difficulties arise when the 
direct definition p.2| ) is used: the convergence of a is often very slow, and it is difficult to tell 
how small the final value of a should be to consider it a = 0. 

A large variety of methods has been proposed to overcome the problem of slowly convergent 
MLCE estimates. Recently, a new algorithm offering excellent convergence, called MEGNO 
(Mean Exponential Growth factor of Nearby Orbits), has been proposed by Cincotta & Sim6| 
( |2000j ). The definition of MEGNO and its mean value is the following ( |Cincotta et al.|2003l ): 

It was shown that if (f){t) is a regular solution with a linear divergence of nearby orbits then 
lim^^oo 0^) {t) = 2, and if (f){t) is a chaotic solution then (Y) (t) ~ {a/2)t, as t ^ (X). More- 
over, when (Y) {t) tends towards a value different from 2, then it indicates that close trajectories 
diverge according to a certain power law. If (f){t) is a periodic solution then (Y) (t) tends to 0. 
The asymptotic behavior of {Y) (t) is given by a uniform formula (Y) (t) ^ at-\-d, where a ~ 
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and d ~ 2 for a quasi-periodic solution, while a cr/2 and d ~ for an irregular and stochastic 
motion. Having Y{t) we can indirectly estimate the MLCE on a finite time interval. The weight 
function s in the definition of MEGNO reduces the contribution of the initial part of the tangent 
vector evolution, when the exponential divergence is too small to be observed relative to other 
linear and nonlinear effects ( [Morbidenil|2Q02| ) . Thus, fitting the straight line to the final part of 
we obtain good estimates of cr from a relatively shorter piece of trajectory than in the direct 
MLCE evaluation. 



3. Modeling the RV data - an overview 

The radial velocity (RV) is still the most efficient technique for detecting extrasolar planets. 
To model the RV signal, the standard formulae by Smart] ( [1949| ) are commonly used. Each planet 
in the system contributes to the reflex motion of the star at time t with: 

Vr{t) = K[cos{u + v{t)) + ecoscj] + Vq, (3.1) 

where K is the semi-amplitude, u is the argument of pericenter, u{t) is the true anomaly (in- 
volving implicit dependence on the orbital period P and time of periastron passage Tp), e is the 
eccentricity, Vb is the velocity offset. We interpret the primary model parameters (i^T, P, e, cj, Tp) 
in terms of the Keplerian elements and minimal masses related to coordinates of Jacobi ( [Lee &| 
|Peale|2QQ3| ) or Poincare HFerraz-MeUo et al.|20Q6 |. 

In our previous work, we tested and tried to optimi ze di fferent tools helpful for exploring 



the multi-parameter space of (x^)^^^ for the model Eq. 3.1 and its generalizations. In the case 
when (x^)^^^ possess many local extrema, we found that good results can be obtained 
with hybrid optimization ( Gozdziewski & Migaszewski'2006'). A single run of the hybrid code 



starts the quasi-global genetic algorithm (GA, [Charbonneau^^ 1 995 j . GA makes it easy to carry 



out a constrained optimization within prescribed parameter bounds or to add a penalty term 
to (x^)^^^- The best fits found with GAs are not very accurate in terms of (x^)^^^' so finally, a 
number of the best fit members of the "population" are refined using a relatively fast local method 
like the simplex of Melder and Nead ( [Press et al.|1992| ). The simplex is a matter of choice, so 
we could use other fast local methods. However, the code using non-gradient methods works 
with minimal requirements for user- supplied information. It is only required to define the model 
function ]the so called fitness function, usually equal to 1/ (x^)^^^] — conveniently, this function 
is the same for the GAs and simplex — and to determine (even very roughly) the bounds of the 
parameters. The repeated runs provide an ensemble of the best-fits that helps us to detect local 
minima of (x^)^^^' ^^en if they are distant in the parameter space. We can also obtain reliable 
approximation to the parameter errors (Bevington & Robinson|2QQ3] ) within the la, 2cr and 3cr 
confidence intervals of x^ at selected 2-dim parameter planes. 

While we prefer the GAs as the quasi-global optimization tool, other efficient and robust al- 
gorithms for identifying and characterizing multiple planetary orbits in precision RV data are 
known. In particular, the Bayesian Kepler periodogram ( Gregory 2007) and Markov chain Monte 
Carlo (MCMC) technique ( [Ford|[2005 ) are proven to be robust tools for calculating the model 



marginal likelihood which is used to compare the probabilities of models with different numbers 
of planets and for investigating the uncertainty of parameters in the orbital solutions. 

Due to the limited time-span of the observations, we often encounter a problem that the data 
only partially cover the longest orbital period. In this situation, it is possible either that (x^)^^^ 
does not have a well defined minimum, or that its shape is very "flat", so the confidence levels 
may cover large ranges of the fit parameters. To illustrate the shape of (x^)^^^ selected 2-dim 
parameter planes, we perform a systematic scanning of the space of initial conditions with the 
fast Levenberg-Marquardt (L-M) algorithm ( Press et al.[p^992) . Usually, for representing such 
scans, we choose the semimajor-axis — eccentricity (a, e) plane of the outermost planet. We fix 
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(a, e) and then search for the best fit, initiating the L-M algorithm with starting points selected 
randomly (but within reasonably wide parameter bounds). The L-M scheme ensures a rapid 
convergence. It is heavily CPU-consuming and may be effectively applied in low-dimensional 
problems. In reward it provides a clear picture of the parameter space. 

In the case of resonant or mutually interacting planets, the problem is even more complex. 
How to interpret the RV measurements in that case often remains an open and difficult question. 
The A/'-planet configurations are parameterized by at least 57V + 1 parameters, even assuming that 
the system is coplanar. The RV signal in terms of Eq. 3. 1 is degenerate — we have no information 
on the inclinations of orbits and the true masses of the companions. Moreover, we do not know 
a priori the number of planets in the system, so the resolved solutions are often not unique. It is 
also well known that the Keplerian (kinematic) model is usually not adequate to properly explain 



the RV variability. Instead, a self-consistent A^-body Newtonian model should be applied ( Rivera 



& Lissauer|2Q01[ Laughlin & Chambers|2Q01j . Basically, the effects of mutual interactions in- 



cluded in the Newtonian model could make it possible to determine or estimate the inclinations 
and masses provided that long enough time-series of precision data are available. Nevertheless, 
the measurements can be affected by many sources of error, like complex systematic instrumen- 
tal effects, short time-series of the observations, irregular sampling due to observing conditions, 
and stellar noise. Little is known on the real statistical characteristics of the error distributions 
and the assumption that these are Gaussian distributions is not necessarily valid (Ferraz-Mello| 
let al.|2QQ5| ), see also ( |Baluev|2QQ7l ). 

All these factors, in particular the unspecified number of planets and undetermined or weakly 
constrained parameters, can (and often do) lead to best-fit solutions representing unrealistic, 
quickly disrupting configurations ( Ferraz-Mello et al. 2005 [ Lee et al. 2006[ Gozdziewski & 



Konacki|^2006). But according to the Copemican principle, the detection of strongly unstable 



systems during two decades of RV observations is not likely, we would rather expect that the dy- 
namical stability should be preserved over a significant part of the parent star life-time counted 
in Gyrs. Stability is therefore a natural requirement of a model solution consistent with observa- 
tions. Many authors take it into account when analyzing the dynamics of the best-fit configura- 
tions with different constraints, e.g., to mention only a few examples in an endless list of refer- 
ences: through long-term integrations, requiring astronomical stability (Lissauer & Rivera*2001 



Laughlin et al.|2002j or Lagrange and/or Hill stability ( Barnes & Greenberg,2007, ), also through 
fast indicators like the diffusion of characteristic frequencies ( Correia et al.|2005 ), fast Lyapunov 
indicators (Dvorak et al.'2005'), maximal eccentricity ( Vogt et al.|2005[ Ford et al. 12005 *), critical 
angles ( ,Ji et al. ^2003^ Lee et al. 2006 ; Beauge et al.^2006^, see also paper by Beauge et al. (2008) 
in this volume), or Te determined over short- time scale (the Systemic project initiated and led by 
Greg Laughlin, www.oklo.org). However, a common approach to search for stable solutions in 
a neighborhood of unstable best-fit configuration by trial and error does not necessarily provide 
stable fits that are simultaneously optimal, in terms of (x^)^^^ or an rms. What is even more 
important, the stability requirement imply non-continuous and complex structure of the space of 
initial conditions that depends on adopted definition of stability. 

Hence, as a general way of modeling the observations we propose to eliminate unstable (for 
instance, strongly chaotic) solutions during the fitting procedure. The idea is very simple: we 
modify the hybrid algorithm by adding a suitable penalty term to the function determining the fit 
quality, e.g., (x^)^^^ or rms for unstable solutions. The penalty term must rely on some signa- 
ture of the system stability. We found that MEGNO is particularly useful for that purpose thanks 
to its rapid convergence and great sensitivity to chaotic motions. We have shown that on many 
examples, this method (called Genetic Algorithm with MEGNO Penalty, GAMP) is very use- 
ful in modeling resonant or close-to-resonant planetary configurations when even small errors 
of orbital phases or other parameters may lead to quick self-destruction of the system. Unfortu- 
nately, the algorithm cannot give a definite answer when we want to resolve the (x^)^^^ shape 
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Figure 1. The levels of reduced (x^)^^^ of the 2-planet Keplerian model of the HD 155358 RV data pub- 
lished in Cochran et al (2007) onto the {Pc, ec) -plane. The left panel is for the whole tested range of 
parameters. The right panel is for the close-up of the most prominent minima around Pc ^ 500 days and 
ec ^ 0.2. The smooth line is for the collision line. Curves labeled with la, 2a and 3a are for the confidence 
levels of the best fit marked with crossed circle. 



in detail or find strictly stable solutions because, in practice, the penalty term can be calculated 
only a over relatively short period of time (due to CPU time requirements). This is particularly 
important for systems affected by long-term resonances, i.e., configurations with three and more 
planets, inclined orbits (Gozdziewski et al.|2007 ). Hence, in an additional step, we need to ex- 
amine the stability of individual best fits selected with GAMP over a period of time related to the 
time- scale of relevant unstable behaviors (resonances). 



4. HD 155358 - a system with two planets 

To illustrate the algorithms and problems discussed above on a new planetary system, we 
consider the RV observations of HD 155358 by [Cochran et al.| ( [2007) . The data published in 
this paper consist of 71 observations and span ^2100 days. The single-planet Keplerian model 
does not fully explain the RV variability, so the discovery team studied a 2-planet Keplerian 
model of the RV. This yields three acceptable fits: a solution which is stable for 100 Myrs with 
Pc ~ 195 days, Pq ~ 526 days and minimizes {xl)^^'^ well as two unstable fits with Pq ~ 500 
and 1500 days, respectively. These fits have very large ec leading to catastrophically unstable 
configurations ( [Cochran et al.|2007) ). 

We try to extend the analysis of the RV data by considering an A^-body model of the observa- 
tions and looking more closely at the phase-space structure of the putative 2-planet system with 
the help of the dynamical tools described above. 

At first, we did two systematic scans of (y?,)^^^ in the (Pc, ec)-plane . The results are shown in 
Fig.[l] Three local minima of (x^)^^^ reported by Cochran et al. ( 2007\ are evident. Two of them 
lie far over the collision line of orbits, defined in terms of semi-axes and eccentricities through 
a5(l + 65) = ac{l — ec). This line marks the zone in which the mutual interactions of relatively 
massive companions can quickly destabilize the configuration. The solution with Pc ~ 300 days 
yields (x?)^^^ similar to that one of the best fit with Pc ~ 500 days. 

Next, we refined the Keplerian fits with A^-body model. The osculating parameters of the 
dominant solution (see Table 1 , fit II) are slightly different from that one of the Kepler model, 
nevertheless it appears also stable. Its neighborhood is illustrated in dynamical maps shown in 
Fig. [2] Curiously, the best-fit configuration it located in the very edge of a stable zone between 5:2 
and 3:1 mean motion resonances (MMRs) of two planets with masses ~ 0.8 mj and ~ 0.5 mj, 
respectively. As we have observed in other cases, the max e and max indicators are in excellent 
correlation with the measure of formal stability (here, log SN). We can also see a very complex 
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Figure 2. Dynamical maps of putative HD 155358 coplanar configuration of two Jovian planets 
( ^ 0.5-0.9 mj). The o sculating elements of the N-body solution at the epoch of the first observation 
in \Cochran et al. | ( 2007\ are given in Table 1 (fit II) and marked with crossed circle. The top-left panel is for 
the Spectral Number. Colors mark the stability regime: black is for regular solutions, yellow is for strongly 
chaotic solutions. The top-right panel is for the maximal amplitude ofapsidal angle /\vj — vjc — vjj^. Pan- 
els in the bottom row are for the max e indicator ( i.e., the maximal eccentricity attained after the integration 
period ~ 30, 000 yr). The most prominent mean motion resonances between the planets are labeled. 



border of the stable zone. The system would be located in dynamically active region of the phase 
space spanned by a few low-order MMRs. The proximity of the best fit configuration to the 
5:2 MMR and moderate eccentricities indicate a dynamical similarity of the HD 155358 system 
to the Solar system. 

To illustrate the hybrid optimization, we performed an independent search for 2-planet Ke- 
plerian fits assuming orbital periods in the range of [100,3000] days and eccentricities in the 
range of [0,0.8]. The ensemble of gathered fits is shown in the top-left panel in Fig[3] To make 
the comparison with the results of systematic scanning more transparent, in this panel we also 
plot contour levels of la, 2cr and 3cr confidence intervals seen in Fig.[T] Note that we plot only 
solutions in the range of Pq G [300, 700] days (compare with the right panel in Fig. [T]). The 
hybrid algorithm also reveals the minima seen in Fig.[T] Besides, it detected one more minimum 
at Pc 300 days and moderate ~ 0.1, at the same depth. This justifies the efficiency and 
robustness of the hybrid code. The algorithm not only detects the best fits solutions but also helps 
to resolve to some extent the shape of (x^)^^^- 

The best fits with Pc ~ 300 days are both very unstable. It does not necessarily mean that in 
their neighborhood some stable solutions do not exist, so it is the case in which the application 
of GAMP can be helpful. Indeed, we can detect two clumps of stable fits (see the top-right panel 
of Fig.[3| in the regime of large ec, both lying far over the collision line. In one of these islands, 
we pick up a rigorously stable solution with (x^)^^^ ~ 1.14 comparable to that one of the best 
fit II (Table 1) yielding only marginally worse rms ~ 6.3 m/s. The evolution of MEGNO and 
osculating elements in this quasi-periodic configuration is shown in Fig. |4] We notice extremely 
large variations of eccentricities up to 0.8. The system would be involved in 5:3 MMR protecting 
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companions from close encounter (see evolution of the critical argument of this resonance in the 
bottom-right panel in Fig.|4]). 

Yet the fit parameters are determined within some error ranges that should be interpreted with 
taking into account the structure of the phase space (see Fig. [2]). To illustrate this problem we 
examined more closely the neighborhood of the best fit II. At this time, we performed two ex- 
periments. In the first search, we applied the hybrid code without stability constraints driven 
by the "usual" A^-body model of the RV. The results are illustrated in the bottom-left panel in 
Fig.|3] The quality of fits within la, 2cr and 3a confidence intervals of the best fit (marked with 
crossed circle; see Table 1, fit II) is color coded with blue, light-blue and gray, respectively. 
The best fits only marginally worse from the best one, are marked in red. Curiously, the plot re- 
veals a subtle structure with three additional local minima of (x^)^^^, in relatively small range of 
ac G [1.1, 1.3] AU. Moreover, these minima are spread over wide range of Cq G [0, 0.7]. Simulta- 
neously, this zone covers many low order resonances, between 5:2 and 3:1 MMR and the (x^)^^^ 
"valley" is crossed by the collision line. Close to this line, the stability could be preserved only 
if the planets are protected from close encounter through an MMR. 

Now, we could examine the stability of every fit that we found but we choose a new search for 
stable solutions in a self-consistent manner with GAMR The results are shown in the bottom- 
right panel in Fig. [3] In this panel, we overplot the stable solutions within la, 2a and 3a con- 
fidence interval of the best stable fit II over all solutions within 3a level which are found in 
the previous search (i.e., without stability constraints). It now is evident that only a part of the 
ixiy^'^ valley can consist of dynamically stable solutions, nevertheless the acceptable fits are 
spread over significant range of Aac ~ 0.2 AU. Apparently, this error is quite small but in fact 
it is large enough to cover a few low-order MMRs. We conclude that the current set of RV data 
cannot fully characterize the system state and new observations are required to constrain the 
elements of the outer planet. 

Finally, both Keplerian and Newtonian 2-planet solutions lead to apparent excess of the resid- 
uals, in particular at the end parts of the RV curve. It may indicate that the 2-planet model does 
not fully explain the RV variability. In particular, the system may involve more than two plan- 
ets. A heuristic argument supporting such a claim may be the proximity of the best fits to the 
collision line. We know similar cases, for instance /i Arae (Gozdzi ewski et al.|2007t PepeetaL 



[2007 ), or HD 37124 ( Vogt et al.|[2Q05] ). To check such hypothesis we looked first for 3-planet 
Keplerian solutions with the hybrid code. The best fit found yields (x^)^^^ ^ 0-87 and signifi- 
cantly better rms ~ 4.6 m/s but is unstable. The GAMP search yields stable configurations with 
quasi-circular orbits of the outermost planets, yielding rms ~ 5.5 m/s. An example fit of this 
type, yielding (x^)^^^ ^ 1-07 and an rms ~ 5.6 m/s, given in terms of osculating element at 
the epoch of the first observation in tuples of (m [mj], a [AU],e, co [deg], M{to) [deg]) is the 
following (0.115, 0.383, 0.025, 281.3, 154.4), (0.770, 0.627, 0.040, 108.9, 173.7), (0.490, 1.187, 
0.000, 359.1, 65.7), for planets 6, c respectively, the offset Vq = 10.21 m/s. We note the small 
mass of the innermost planet. Its RV signal is at the level of noise, so additional observations 
would be required to confirm of withdrawn such a model. 



5. Trojan planets in the Gruis system? 

[Laughlin & Chambers | ( [2002] ) predict that the reflex signal of a single planet in a quasi-circular 
orbit may be also interpreted by two Jovian Trojan planets, i.e., two objects sharing similar orbits 
(involved in 1:1 MMR). That possibility is intriguing because stable Trojan companions to the 
stars may be quite common. It can be indicated by a number of stable Trojan configurations in the 
Solar system. Some argue that they can be a frequent by-product of planet formation and and/or 
dynamical evolution (Laughlin & Chambers 2002 ). However, the genesis of Trojan planets is not 
quite clear because on the contrary, there is some evidence that formation of such bodies could 
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Figure 3. The top-left panel is for the (xlV^'^ of the 2-planet Keplerian solutions to the HD 155358 RV 
data published in Cochran et al. (2007) derived with the hybrid algorithm. The fit parameters are projected 
onto the (Pc, ec)-plane. Their quality is color coded: dark blue is for la-, light-blue is for 2a- and grey is 
for 3a-confidence levels of the best fit marked with crossed circle. Contours are for the confidence levels 
obtained in the systematic scan (see the right panel in Fig. 1). The top-right panel is for stable solutions 
in the range Pc ^ 300 days derived with GAMP driven by the N-body model. The stable best-fit solution 
in this area is marked with crossed circle, and its elements are given in Table 1, fit I. The bottom-left panel 
is for the ensemble of fits gathered with the hybrid algorithm driven by the N-body model of the RV data 
around the dominant minimum of (xlV^'^- The bottom-right panel illustrates stable solutions gathered in 
the GAMP search. The set of solutions within 3a of the best stable fit (elements given in Table 1, fit II) but 
obtained without stability constraints, as shown in the bottom-left panel, are again shown as a gray-filled 
contour. Approximate positions of the most prominent MMRs are labeled. 



be difficult ( Beauge et al.||2007| ). Still, many authors expect that Trojan planets can exist [see, 
for instance, the work of Dvorak et al. in this volume and references therein, also ( [Ford & Gaudi| 
[2006) ]. 

Recently, we found a similar kind of ambiguity of the RV models concerning 2:1 MMR con- 
figurations. At present, we know five extrasolar systems presumably involved in 2:1 MMR, i.e., 
Gliese 876 (Marcy et al 2001) , HD 82943 (Mayor et al..2003), HD 128311 (Vogt et al. 2005])^ 



HD 73526 ( [Tinney et al.|2006| ), and /i Arae ( [Jones et al. 12003} [Gozdziewski et al.|2007t|Pepe et al. 



2007| ). However, the 2:1 MMR model of the radial velocity observations can also be non-unique. 



The periodogram of the 2: 1 MMR RV signal is very similar to that one of the 1 : 1 MMR. Indeed, 
the RV variability of HD 128311 and HD 82943 can be explained by highly inclined systems 
in 1:1 MMR (Gozdziewski & Konacki|2006| ). We also found that the RV of HD 73526 can be 
modeled with two highly inclined Jovian Trojans. The modeling of the 1 : 1 MMR is a challenging 
problem because Jupiter-like planets sharing eccentric orbits with similar semi-major axes inter- 
act heavily and the collisional configurations are generic. Hence, stability constraints are critical 
in the search for optimal and stable configurations. This seems to be one of the best applications 
of GAMP like algorithms. 

Among a few cases we analyzed so far, the Gruis appears to be a particularly interesting ex- 
ample of the possible "Jupiter on circular orbit"-"two Trojans" ambiguity. A Jovian companion 
to the GO dwarf Gruis in a wide and almost circular orbit has been announced in the work by 
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Figure 4. The temporal evolution of orbital elements in the 2-planet, Newtonian solution to the RV data 
ofHD 155358 related to highly eccentric orbits (see the top-left panel in F/g. The osculating elements 
are given in Table 1 (fit I). The top-left panel is for the semi-major axes, the top-right panel is for the 
eccentricities. The MEGNO is plotted in the bottom-left panel The bottom-right panel is for the critical 
angle of the 5:3 MMR. 



Jones et al. ( 2003|). In our analysis, we use updated RV data comprising of 59 precision measure- 
ments ( [Butler et al.|2Q06| ). The best-fit single planet model to these data yields ^ 1300 days, 
and ~ 0.1. We re-analyse the data to look for possible Trojan planet solutions. Curiously, we 
quite easily found many stable, coplanar configurations involved in 1:1 MMR (Fig. [5]) yielding 
similar or slightly better fit quality (rms G [5, 6] m/s). The reflex signal of the 1:1 MMR (Fig.|6]) 
can hardly be distinguished from that of a single-planet system. The osculating elements of the 
Trojans are given in Table 1 (fit III). In this case, both planets would move on quasi-circular 
orbits and these would be coplanar. The dynamical maps shown in Fig. [5] (also accompanying 
other best-fits solutions with acceptable quality which we found in the search) illustrate the ex- 
treme variability of the 1:1 MMR islands. The map for the best-fit with elements in Table 1 is 
shown in the top-left panel of Fig. [5] Another peculiar solution is illustrated in the top-right map 
in Fig. |5] The initial eccentricities are moderate, and the system would be found in extremely 
large island of stable motions. It spans whole range of e\^. The possibility of existence of such 
extended stable zones may strength the hypothesis of stable extrasolar Trojans. 

The ambiguity of the RV fits implies interesting issues concerning the models of creation and 
stability of Earth-like planets interior to the orbits of the putative Jovian Trojans. In the Gruis, 
the space interior to the Jovian planet is "empty" as no smaller planets have been yet detected. 
So we can try to predict in which regions of the habitable zone (~ 1 AU) smaller planet could 
survive. For this purpose we computed dynamical maps for putative Earth-like masses with initial 
conditions varied in the (ao, eo) plane, and initial orbital angles set to 0°. We considered two 
dynamical environments: the one with the best-fit Jovian companion in close to circular orbit and 
the second one with Trojans in quasi-circular orbits (their elements are given in Table 1, fit III). 
The results are shown in Fig|7] For the first configuration, we detect an extended zone of stable 
motions. Additional experiments regarding creation of Earth-like planets through coagulation 
of Mars and Moon- size protopolanets (see Raymond, 2008 in this volume) performed with the 



Mercury code (Chambers||1999| assures us that such planets emerge easily in that zone. In the 



case of a configuration with Trojans, the stable zone shrinks significantly. Moreover, the creation 
of Earth-like planets is much more difficult. We found that they could form only in the zones of 
relatively stable motions, up to ~ 0.9 AU and in the "gap" between the 4:1 MMR and the border 
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Figure 5. Dynamical maps of putative Gruis coplanar, edge-on configuration of two Jovian planets 
0.5-0.9 mj) involved in low-order resonances. The best-fits yield an rms ~ 5-6m/sand (xlY^^ ^ 1- 
Their quality is similar to that of the single-planet solution (an rms about of 6 m/s). 
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Figure 6. The RV data ofr^ Gruis and the synthetic signals. The red curve is for the the best, single-planet 
Keplerian fit yielding (x^)^^^ ~ 1 and an rms ^ 6.1 m/s. The blue curve (darker one) is for coplanar, 
edge-on configuration (an rms ^ 5.9 m/s) involved in lb:lc MMR (see the top-left panel in Fig. |5[/or the 
dynamical map). 

of global instability. Yet in that case, the simulations are very difficult to carry out due to frequent 
close encounters between planetesimals and the Jovian planets. 



6. Conclusions 

In this work we consider some problems related to modeling observations of stars hosting 
multi-planet systems. It is well known that the phase space of such system has a non-continuous 
and complex structure with respect to any stability criterion. Hence, when searching for initial 
conditions, one has to take into account the dynamical character of putative planetary configu- 
rations. Due to narrow observational windows, significant measurement errors, stellar jitter and 
other uncertainties, the formal best-fits may appear very unstable. Searching for stable solutions 
in their neighborhood of the phase space by trial and error, we should not expect that the results 
could be statistically optimal. Thus, an intuitively natural approach is to eliminate unstable con- 
figurations during the fitting process, through penalizing unstable solutions with a suitably large 
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Figure 7. Dynamical maps for Earth-mass planets in the coplanar Gruis system. The left column is for 
the stability map for the best fit configuration with one Jovian planet in quasi- circular orbit (a ^ 2.5 AU), 
the right column is for the systems with Trojans (Table 1, fit III). Panels in the top row are for the stability in- 
dicator log SN. Panels in the bottom row are for the max e indicator ( the integration period ^ 50, 000 yr). 
Some MMRs between the Earth-like planet and the Jovian companions are labeled. 



HD 155358 (fit I) HD 155358 (fit II) Gruis (fit III) 

Parameter planet b planet c planet b planet c planet b planet c 



m sin i [mj] 


0.827 


0.490 


0.863 


0.497 


0.401 




0.923 


a [AU] 


0.623 


0.875 


0.628 


1.212 


2.471 




2.565 


e 


0.121 


0.743 


0.128 


0.198 


0.027 




0.053 


cj [deg] 


131.6 


86.11 


161.9 


272.0 


99.9 




163.8 


M{to) [deg] 


106.32 


313.7 


130.3 


198.8 


3.5 




3.6 


1.14 




1.08 






1.12 




Oj [m/s] 


5 




5 






4 




rms [m s~^] 


6.31 




5.98 






5.96 




Vo [ms-^] 


12.69 




10.69 






-0.04 




M. [Mo] 


0.87 




0.87 






1.25 





Table 1. The best-fit astro-centric, osculating Keplerian elements of stable, coplanar and edge-on planetary 
configurations at the epoch of the respective first observation. Original errors of the data are rescaled by 
adding the "jitter" aj in quadrature. 

value of the (x?)^^^ function, or of another measure of the fit quality. In that way, the stability 
plays a role of an additional, implicit observable. That method is suitable for multi-body systems 
with Jovian planets presumably involved in low-order mean motion resonances (IMlVERs). In par- 
ticular, we considere d two new examples in which the model of the RV may be non-unique. The 



mples in 

RV of HD 155358 by|Cochran et aLlhoOV*) permit a few local minima of {xlV^^ related to dif- 



ferent orbital configurations. We also found an example illustrating the ambiguity of Keplerian, 
close to circular single-planet solutions. The RV data of Gruis could be equally well modeled 
with coplanar configurations of Jovian planets involved in 1:1 IMlVERs. 

IMoreover, our fitting method, used mainly for RV data, is quite general and may be applied to 
other types of observations as well. As the stability criterion, one can use the maximal Lyapunov 
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exponent, the most stringent and formal characteristic of stable/unstable motion. Other suitable 
indicators like the maximal eccentricity, the spectral number, or the diffusion rate of fundamental 
frequencies may also be applied. These fast indicators help us to search for and find long-term 
stable solutions, but also make it possible to efficiently explore and to visualize the sophisti- 
cated and varying structure of the phase space. We can see the planetary system in its dynamical 
environment. 

Many multi-planet systems are found on the edge of long-term dynamical stability. It is not 
clear yet whether this is a general property of multi-planet systems, the outcome of poor statistics 
or just the consequence of a bad choice of the RV model. Large eccentricities in multi-planet 
systems may "hide" other, unknown planets. Yet in that case, the dynamical modeling of the RV 
with stability constraints provides valuable information on the dynamical structure of the putative 
planetary configurations. Finding the best fits on the very edge of stable zones may provide good 
hints and motivation to look for alternate models of the RV. 
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